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ABSTRACT 

For massive black hole binaries produced in galactic mergers, we examine the 
possibility of inspiral induced by interaction with field stars. We model the evolution 
of such binaries for a range of galaxy core and binary parameters, using numerical 
results from the literature to compute the binary's energy and angular momentum 
loss rates due to stellar encounters and including the effect of back-action on the 
field stars. We find that only a small fraction of binary systems can merge within a 
Hubble time via unassisted stellar dynamics. External perturbations may, however, 
cause efficient inspiral. Averaging over a population of central black holes and galaxy 
mergers, we compute the expected background of gravitational radiation with periods 
Pu, ~ 1 — lOy. Comparison with sensitivities from millisecond pulsar timing suggests 
that the strongest sources may be detectable with modest improvements to present 
experiments. 



Subject headings: black holes — galaxies:interactions — pulsars:timing — radiation 
mechanisms : gravitational 



1. Introduction 



The possibility that many galaxies contain a central black hole of 10 Mq or more, coupled 
with the ideas that elliptical galaxies probably formed from mergers (c./. Barnes & Hernquist, 
1992) and that all galaxies may have formed from multiple mergers of subunits, has led naturally 
to consideration of the fate of two or more massive black holes in a merger remnant. Simulations 
suggest that the cores of victim galaxies merge via violent relaxation in a few crossing times 
to form the remnant core (Barnes & Hernquist); with some possible exceptions (Governato et 
al. 1994), black holes in the victims then end up in the core of an elliptical or a spiral bulge. 
The holes sink towards the center of the stellar distribution on the dynamical friction time scale 
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of ~ 10^ years until they become bound to one another (Begelman, Blandford k, Rees 1980, 
hereafter BBR), forming a massive black hole binary (BHB). If the binary loses enough energy and 
angular momentum to the field stars, which is by no means certain, it will enter a regime where 
gravitational radiation alone can bring about inspiral and coalescence within the Hubble time. 

During BHB inspiral, the gravity waves generated sweep through a range of frequencies. In 
the final merger of two 10^ black holes, a 'chirp' of gravitational radiation with periods as 
small as Pqr ~ lO^s will be emitted. It has been suggested that this radiation might be detected 
through the Doppler tracking of spacecraft and by using space-based interferometers (Thorne & 
Braginsky 1976, Thorne 1992, Fukushige et al. 1992a). Unfortunately, the binary spends only few 
wave periods in the final chirp, so if inspiral is rare such events will be difficult to observe. In 
contrast, at binary frequencies ~ 0.01 — 1 /^-Hertz, the BHB persists for many orbits. As first 
pointed out by Detweiler (1979), the strain amplitudes of the ultra-low frequency (ULF) gravity 
waves generated by sufficiently massive BHBs may be detected out to cosmological distances as 
perturbations in the pulse arrival times of quiet pulsars (Romani &; Taylor 1983). With a new 
generation of timeable millisecond pulsars (Kaspi, Taylor k, Ryba 1994, hereafter KTR; Camillo, 
Foster & Wolszczan 1994) allowing remarkable improvements in sensitivity, it seems appropriate 
to assess the probability of massive BHB detection. 

In this paper we examine the fate of BHBs in a population of merger remnants, estimating 
the frequency of events and using these results to outline the required pulsar timing sensitivity 
needed to effect a detection or constrain the merger hypothesis. First we consider inspiral as the 
average result of repeated interactions with the surrounding stars (Mikkola k. Valtonen, 1992, 
hereafter MV; Hills, 1993). The binary loses energy, shrinks and moves faster (hardens), with 
moderate increase in eccentricity. Numerical simulation (Makino et al. 1993, hereafter MFOE) 
seems to confirm this evolution in general, but warns of the importance of back-action of the BHB 
on the surrounding core, neglected in previous analytical investigations. Since the characteristics 
of the host core are crucial to the success of the inspiral, in section 2 we develop a numerical 
integration of the evolution of binding energy and angular momentum (hence eccentricity) during 
a merger event, which can be applied to a range of black hole masses and core properties. This 
approximate trajectory includes back-action on the core via loss cone depletion (BBR) and 
heating. In Section 2.3, we estimate fsi, the fraction of mergers achieving a successful inspiral due 
to stellar scattering, by computing our evolutionary model for a range of core parameters and BHB 
properties. Additional mechanisms that may drive inspiral are also mentioned in Section 2.3. In 
section 3 we collect estimates of the number density of massive black hole cores and galaxy merger 
rates to estimate the number of merging sources visible as a function of gravity wave period and 
strain amplitude. In section 4 we compare these rates to present and anticipated pulsar timing 
sensitivities and conclude by commenting on the prospects for further constraining inspiral events. 
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2. Simulation of Inspiral 



2.1. Binary Evolution Formulae 



It was suggested by Fukushige et al. (1992b) that the evolution of a BHB can be computed 
from the classic Chandrasekhar formula (Chandrasekhar, 1943) for dynamical friction (DF) on 
a single object in a homogeneous background. Applying an approximate formula to both holes, 
they find a large growth in the eccentricity; this was our motivation to trace energy and angular 
momentum loss of a BHB. We use the results of MV, who track the net effect of an ensemble 
of three-body interactions between an equal-mass binary and an intruder from a background of 
light objects. They cast their results in the form {a/a?) = {Ra/v^)i:Gp, where a is the binary 
semi-major axis, p the density of background stars, and v^o the pre-encounter intruder speed. 
They find the numerical factor Ra to be well fit by Ra = — 6.5/[l -|- 2.44(tioo/H^)^], independent of 
binary eccentricity. Here W is the r.m.s. relative speed of the binary elements, i.e. the relative 
speed of a circular binary of equal energy. 

We assume the galactic core to have a Plummer model velocity distribution, given by 



Since the one-dimensional dispersion a depends on position in a Plummer model, we use the 
central value ctq. Using the fitting function, we calculate (as a function of do) {o'oRa/'Voo)aoy the 
average over Plummer distributed Voo- This gives 



circular binary due to DF, applying the Chandrasekhar formula (c.f. Binney &: Tremaine, 1987) 
independently to each hole. In the W < ao limit, the MV formula agrees by construction with the 
analytic calculation of Gould (1991); it thus differs from the DF result in the form of the Coulomb 
logarithm A (see Gould for a discussion). For the DF calculations of Fig. 1, we take the maximum 
impact parameter bmax contained in A to be a function of the binary period, as discussed in the 
next section. In the W ao regime which proves most important in our calculation, we see in the 
figure that DF does not approximate well the more physically grounded and asymptotically sound 
three-body results. 

For the eccentricity evolution, MV give results in the form 




(1) 



for the evolution of the binary energy. 



E = 




(2) 
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Xn+l/2 _ X 



where X = 1 — 
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For each X and W we again average over with distribution (^), and convert (^) to a time 
derivative, giving 



^n+l/2 _ 



(4) 



In Fig. 1 we also compare the amount of eccentricity growth of the binary as a function of energy 
for DF and MV. Again, DF does not seem to be a suitable approximation. Since MV treat only 
the equal-mass binary, we impose a correction factor to force agreement with Gould's analytic 
mass ratio scaling for Mi ^ M2. We use Mi > M2 throughout. 
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Fig. 1. — Comparison of binary evolution formulae given by dynamical friction and by Mikkola and 
Valtonen (1992). The left scale applies to the solid lines giving E, and the one value in this range 
from Hills (1983); the right scale is for the eccentricity growth rate at eccentricity of 0.9 (dotted 
lines). To calculate bmax and A for the DF formulae we use Mi = M2 = I.OM9, Tc = 200pc and 
a = SOOkms"^ 



In following the orbit evolution, we also include the energy and angular momentum losses to 
gravitational radiation (GR): (Peters 1964) 



/ ( 



dE_\ _ 32 G^M{Mi{Mi + M2) A 73^2 



\ dt 



e2)7/2 
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dL\ _ 32 G'^/^M^ Ml {Ml + IVh)^^^ f 7 2\ 



(5) 



(6) 



Once these terms dominate, the BHB will circularize. When the time scale for binary evolution 
by GR alone is less than 10^'^ years, we call this the GR regime. 
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2.2. Properties of the Core 



We use a Plummer model core of stars with mean mass m^, core radius and isotropic 
central velocity dispersion (|^). The core stars have total potential energy W and kinetic energy 
K, satisfying the virial relation 2K + W = 0. With 

W{N r)- GjNsmsf GNsmsjMi + M2) 
32 Tc 

we use K{Ns,ao) = {9'7T/32)NsmsaQ and the virial relation to find 

^^,^_6agr, 32 Mi + M2 ^ 

The core mass is Mc = Ngms- 

Only a subset of the core stars can effectively exchange energy with the binary. According to 
Heggie (1975), these are the stars whose interaction time with the binary is less than its orbital 
time scale. This criterion defines a maximum impact parameter bmax (a function of the binary 
period) which we also used in the DF calculation for Fig. 1. At any distance from the center of 
the core, the condition on impact parameter selects stars whose velocities are within a cone aimed 
towards or away from the center, the so-called loss cone (Frank & Rees 1976). For our core, the 
fraction of stars with impact parameter at most bmax is found by direct integration to be 



^avail _ ^max /q\ 
s "max I ' c 

Since scattered stars enter a solid angle at a rate proportional to the solid angle size, stars 
can enter this loss cone no faster than the rate N avail /treiax- Thus tre/aa; is also the loss cone 
repopulation time, 




tic = trelax = 7^ 3.7 X lO^^yr — (10) 



f f \ / SOOkms"^ 



VlOiOy V200pc/ \^ do 



(Binney & Tremaine 1987). 



Once encounters with the binary deplete the loss cone (BBR), the binary evolution rate will 
be limited by the rate at which stars enter the loss cone, and the amount of energy each star 
can take from the binary. In the W > limit, MV find the change in binary energy with each 
stellar encounter to be {5E) ~ 2Ems/{Mi + M2). This is in good agreement with the results at 
several mass ratios of Roos (1988) and Hills (1983). Writing this as {6E) = {M2 / Mi)msV2c, where 
the circular velocity V2c of the lighter binary component is roughly the closest approach velocity 
of a star, we pass to the W < regime by using (Tq for V2c- The energy the loss cone can take 
away without repopulation is thus Eavail = Navaili^E), and from (|lO|), energy enters the loss cone 
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at a rate {dE/dt)in = E avail /treiax- If each star which interacts with the binary is scattered into 
a random orbit, a fraction Navaii/^s of stars wih be returned to the loss cone. For a loss cone 
population in steady state, then, the limiting binary evolution rate is 

fdEbin\ ^ Ng ('^^\ = ^"""^ {5E) (11) 

\ dt J -^s ^avail \ dt J ^relax 

This rate is generally much smaller than that due to unrestrained dynamical friction. We 
approximate the situation by using the average core density in (Q) and (^) until Eavaii drops below 
the binding energy of the binary. From this point on the binary will affect the loss cone stellar 
density significantly, and we impose the steady-state evolution rate. Invariably this 'loss cone 
catastrophe' results in a binary evolution time scale Ehin/Ehin much longer than the Hubble time. 
Similar but less specific treatments have been presented by BBR and Roos (1981). 

We also treat eviction of stars from the core, found to be significant in the N-body work of 
Makino et al. (1993, hereafter MFOE). We assume all stars which interact with the binary after 
V2c is larger than the Plummer central escape velocity v^sc = V^o'o receive enough recoil velocity 
to be kicked out of the core. We calculate the resultant change in the core energies assuming rapid 
revirialization, readjusting Tc and ctq. The binary also heats the core before becoming hard, but 
this effect is small. 

At the beginning of the simulation, the binary's semi-major axis is the binding radius r;,, 
defined such that the mass of stars within is less than Mi. Its center of mass is at the core 
center, and it has some initial eccentricity Cj. We assume the core is not rotating after the merger, 
since N-body merger simulations show very small angular momentum in the remnant bulge 
(Barnes 1992), even when progenitor bulges rotate (Hernquist 1993). A simulation is terminated 
unsuccessfully if Ej^n/Ebin exceeds 10^*^ years, or successfuly if the GR regime is reached. To 
estimate the parameter space available for successful inspiral, we find log(l — ef) for the marginally 
successful trajectory to one part in a thousand by using ten bisection steps on the initial value. 
Some sample trajectories are shown in Figure 2. If we ignore back-action and loss cone depletion 
the results are identical to those of MV. 

Our simulations confirm the basic behavior described by MFOE. For a range of initial 
conditions, interaction with background stars can cause some increase in the binary eccentricity. 
If we choose core parameters and black hole masses to match those of MFOE, we find loss cone 
depletion to occur at binary radii similar to those at which they terminate their simulation, where 
they too note a slowdown in binary evolution which they attribute to a possible loss cone effect. 



2.3. Success Rate and Alternative Merger Mechanisms 

For all pairs of binary member masses drawn from the set logM/M© = {10, 9.5,. ..,7.0, 6.5} 
we calculated the average initial critical eccentricity required for successful inspiral in 80 galactic 
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Fig. 2. — Trajectories for two BHB/cluster combinations and several Cj. Slight eccentricity growth 
and termination at the loss cone limit are visible, as is recircularization in the GR regime. Binding 
radius r^, and the radius at which V2c > co are indicated. 

cores, drawn from a smoothed distribution based on the data set of Lauer (1985). We chose the 
distribution by noting that the binned data show N{M/10^Mq) = exp(— M9/22) (presumably 
the tail of a Schechter function); matching this and the N{rc) of the data set gives our core 
distribution in the rc — ctq plane. In drawing from the distribution we required that the stellar 
mass of the host core be at least twice the total mass of the holes. 

To estimate the fraction of successful inspirals fsi for any given binary mass, we assume 
stochastically distributed initial eccentricities {P{e) = 2e), so the inspiral success probability 
is P{e > Ccrit) = 1 — ^Irif This choice of probability may be somewhat too low since N-body 
simulations of collisions of bulge-disc-halo systems (Barnes 1992) show radial final infall of the 
bulges; high initial binary eccentricity seems a likely consequence. Some critical eccentricities 
for BHB (log log -^) which will later prove most relevant for pulsar detection are: 
(9.5, 9) ^ .985, (9, 9) ^ .987, and (9, 8) .990. These yield fsi of 2.9%, 2.6% and 2.0% 
respectively. Over the whole range of binary masses, the fsi are ~ 1% — 3%. From the dotted 
curves in Fig. 1 we see that the dynamical friction formulae yield much more eccentricity growth, 
hence earlier arrival in the GR regime as per equations (^) and dH), and larger success fraction. 
We ran calculations with all the above features but using DF, finding fsi ranging from 1% to 
50%. In the case of DF, it is possible to account for rotation of the stellar core, which was found 
to be an important impediment to eccentricity growth in the simulations of MFOE. By giving all 
angular momentum lost by the binary to the core, and applying the DF formulae to the motion of 
the binary relative to the rotating background, we find much less eccentricity growth, resulting in 
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fsi between 1% and 16%. If even the modest eccentricity growth in the MV picture, as seen in 
Figure 2, is suppressed by spinup of the core, the MV fsi could be driven slightly lower. 

Clearly stellar action on an isolated BHB in a galaxy core is of limited efficacy. However, the 
absence of many double active galactic nuclei (AGN's) and of orbital variations in the broad line 
region, as well as short binary periods inferred from from AGN jet precession (Roos 1988) suggest 
that inspiral beyond the loss cone limit does occur in many cases. Other means of binary inspiral 
have been proposed: gas may be ejected during nuclear activity, taking with it energy and angular 
momentum acquired from the binary; or it may be accreted by the larger hole, causing orbital 
contraction as Mir 2 remains constant to conserve angular momentum (BBR). 

Roos (1988) has also pointed out that if merger events happen several times during the 
lifetime of a galaxy and if, in particular, mergers are related to the onset of accretion and nuclear 
activity in AGNs, then external perturbations due to incoming galactic masses should accelerate 
stellar scattering inspiral. ThoTigh a binary which fails to merge is left evolving on the relaxation 
(and loss cone filling) time scale, Roos estimates that by the time a new merging galaxy is within 
a few core radii of the remnant pair, the loss cone can be refilled quickly and repeatedly, driving 
the original binary to GR inspiral. His calculation of the flux of stars entering the loss cone in 
this tidal repopulation scheme suggests that the galactic orbit decays no faster than that of the 
binary. If the incoming galaxy also has a black hole, it seems likely to form a binary with the 
post-coalescence hole. 

Accordingly, in the remainder of this paper we follow two hypotheses: that inspiral occurs 
with the low probabilities allowed before loss cone depletion by MV's stellar encounter dissipation, 
and that repopulation ensures that inspiral always occurs. The latter picture assumes a galaxy 
formation model involving repeated mergers; these two cases certainly bracket the actual situation. 



3. Merger Rates and Populations of Nuclear Black Holes 

We wish to compute the rate of mergers observed from some redshift z. In an Einstein - 
de Sitter universe (A = 0, = 1) with the metric ds^ = dt'^ - a{t)'^[dx^ + x^dJl^], let Fmiz) 
be the number of merger events at redshift z in the history of today's bright (L*) galaxies, per 
dimensionless comoving volume per unit redshift. With this definition, / Fjn{z)dz = NugaiaQ, 
where Ugai is the current number density of such galaxies and oq the current scale factor; N can 
range from ~ 1/3, if only single mergers occur to form elliptical galaxies, to ~ 10 in scenarios 
where all galaxies are formed from repeated merging of building blocks. At redshift z we observe 
a dimensionless co-moving area Anx^iz), with x{z) = ^^^^^ (^l — ^ ■ During a period at 
redshift z, the comoving volume observed is Vy^ = 47rx^(c/a)r2. Thus the number of mergers per 
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unit observers' time is 

K^)d. = ^^f^ = levr f^)' 4 (l - ^) V^(z)d.. (12) 

Note that oq vanishes in any physical rate because of the normahzation of Fm{z). 

The form of Fm{z) is still quite uncertain. Burkey et al. (1994) find that the population of 
close galaxy pairs which seem certain to merge varies as ~ (1 + z)^'^^'^, in the interval < z < 0.6. 
Colm et al. (1994) find that a galaxy density ~ (1 + z)'^'^ best fits the total galaxy count data in 
a model which accounts for photometric evolution. The 3.5 power law gives Fm{z) ~ (1 + z)'^'^. 
We normalize Fm{z) to the merger rate implied by the pair counts in Burkey et al. Applying this 
rate in the interval < z < 1 implies 40galaxies will have suffered a merging event; alternatively if 
this rate continues back to the epoch of high quasar activity, z ~ 3, then each bright galaxy will 
have experienced roughly 5 merger events. We find: 

iy{z) dz = 7.6 X 10-2 yr-i/i5o' f TT^rS^— 3 ) (1 + ^)'-' fl " -7=] ' d^, (13) 
where Hq = 50/i5okms^"'^/Mpc. 



io-3Mpc-3y ^ ' \ ^/r+ 



As an alternative less dominated by early merging, we take merging rate per comoving volume 
constant in time, so that Fm{z) ~ (1 + z)-^/^. If we assume that merging began at some Zm and 
continues to the present, resulting in N mergers per bright galaxy of number density Ugai, we find 

.(.)dz = 0.55yr j (ipi^ j ( + ,^)3/2 _ 1 j ( (1 + ,)5/2 j ^z (14) 

If accretion onto massive central black holes is the source of AGN luminosity, the population 
of remnant holes can be estimated from models of AGN evolution (Cavaliere and Padovani, 1988; 
Small and Blandford, 1992). Recent HST detections of kinematic evidence for massive compact 
objects in nearby galaxy cores support this scenario. We adopt here Small and Blandford's more 
conservative model lA, with a flat local Seyfert luminosity function. This gives the black hole 
number density spectrum 

7.0 X lO^^Mg-^dMg Mpc-3 logMg < -1.6 



Nii{Mg)dMg = { 1.2 X lO-^Mg-^-^dMg Mpc'^ -1.6 < log Mg < 0.45 (15) 

^9" 



1.1 X lO^'^Mg-^-^dMg Mpc-3 0.45 < logMg 



Comparing the above to a Schecter luminosity function of bright galaxies 
N{L)dL = 2.5 X 10-3^|o(^/^*)~^'^e-(^/-^*)d(L/L*)Mpc-3 (de Lapparent et al. 1989), we 
found that if all massive black holes are in bright galaxies, the integral of equation (^) implies 
that 21% of L > L^< galaxies contain a black hole with M > IO^'^Mq today. Gonservatively, we 
hold the co-moving number density of holes fixed: if there are Nm subunits at redshift z that will 



- 10 - 



become a bright galaxy, then the number of BHs per subunit is We note that the high 

inspiral rates (~ 1/3 — 10yr~^) estimated by Fukushige, et al. (1992a) are based on excessively 
optimistic assmnption of an M > 10^ Mq black hole in each of ~ 10 galaxy subcomponents forming 
an elliptical. In their picture there are 10 gravity wave chirps emitted for each bright elliptical 
seen today. As we shall see, such event rates are not consistent with bounds from pulsar timing. 



4. Gravitational Wave Amplitude of BHB Inspiral 

For a circular binary of orbital period Pb, reduced mass and total mass M, Thorne (1987) 
gives a characteristic strain amplitude averaged over direction and polarization of 

he = 8(2/15)^/V(2vrM/Pb)2/3/r (16) 

for the emitted gravitational wave from a circular (e = 0) binary. In units convenient to the BHB 
problem this is 

= 1.8 X 10-15MiM2(Mi + M2)-i/V^;j5'/2(g^^)/^]p-2/3 

where masses are in units of lO^MQ, distances are in kiloparsecs, g{e,n) gives the fraction of the 
wave power in the n^^ harmonic as a combination of Bessel functions (Peters & Mathews 1963) 
and the observed gravitational wave period in years is = Pb{l + z) jn. Only the n = 2 harmonic 
is non-zero for a circular orbit. 

Since the gravity wave flux scales as h?uj~'^ , falling off with the luminosity distance as d^^, 
one can write the gravity wave amplitude from a source at z emitting waves that ha,v6 period. Pyj 
at redshift in an Einstein-deSitter universe: 

/i_i5(z) = 0.29MiM2(Mi + M2)-i/=^/i5o[5'/'(e, n)/n]p-^/ \ ^K\/2^ (1^) 

[1 - (1 + z) ^/^J 

in dimensionless units of 10~^^. The characteristic lifetime for the source at this period is 

TGR = 1.3 X m^M^^M^\Mi+M2f'''f{e)-\nP^/2{l + z)f^y (19) 
where /(e) = (1 + i + § e^)/(l - 6^)7/2, as in Eq. (|). 

We wish to estimate, for a set of assumptions about the population, the expected 
number of BHBs detectable at a given strain sensitivity /i_i5. Inverting equation ( p^ gives 
z{h; Ml, M2, Pw,e)- We can combine this with the merger frequency rate, two independent draws 
from the black hole mass function ([T5|), the successful inspiral fraction for these hole masses 
(averaged over cores that might contain such holes), and the lifetime of the resulting binary at the 
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indicated orbital period to get the number of visible gravity wave sources in a given amplitude 
range: 



diV(fe,Mi,M2,P^,e) 
dMidMa 



iy{z)dz {fsi{Mi, M2)}tgr{Mi, M2, Pw, e, z) — . (20) 



Integration of this equation over Afi, M2 and a range of h gives the source population estimates 
shown in Figure 3a for waves of period Pyj = lOy. In practice, we find that even for inspiral driven 
by stellar encounters, the characteristic eccentricity of the binary decreases to ^ 0.3 by the time 
it reaches the range observable with pulsar timing, so restriction to circular orbits is reasonable 
for estimating fluxes. Depending on assumptions, we see that gravity waves should be detectable 
once a strain sensitivity of ~ 10~^^ is reached. The properties of the detected BHBs can also be 
estimated by weighted sums over this population; for the case of fsi = 1, and mergers proceeding 



as in equation (IS), average characteristics of the detected binaries are shown in Figure 3b. It is 
also worth noting that with these parameters, we expect ~ 1 chirp per lO^yr with an amplitude 
greater than 10~^^ at periods Pw ~ lO^s, substantially lower than the rates quoted by Fukushige 
et al. (1992a). 
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Fig. 3. — Left - Source intensity distributions for several merger models. The filled dots give the 
distribution for tidally repopulated merging and open circles show numbers for evolution slowed by 
loss cone depletion, with the merger rate (13). Crosses show the population from mergers constant 
in t. Right ~ weighted observed source properties for fsi = 1 mergers following rate (13). The 
binary mass ratio is q, and the GR inspiral time is r. 



5. Conclusions and Observational Prospects 
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Millisecond pulsar timing has been used to place strong bounds on the energy density in a 
stochastic background of ultra-low frequency gravitational waves (KTR and references therein). 
For example, these studies have placed an energy density limit Qg < 2 x and have 

helped to rule out various exotic cosmologies, such as structure formation seeded by cosmic string 
loops. With a number of high-quality pulsars now being timed, it is worth considering whether 
astrophysical sources, such as BHBs, are within reach. In Figure 4 we show a simulated gravity 
wave spectrum computed from the amplitude number distributions (as in Figure 3a) for a range 
of wave periods. Two models are shown, the first with inspiral mediated by tidal repopulation 
{fsi = 1), the second for evolution limited by loss cone depletion. Amplitude distributions were 
computed at each frequency; then the power was integrated over the low h sources and Monte 
Carlo sampled from the rare bright sources. The corresponding amplitude is given in Figure 4, 
which thus shows the expected BHB ultra-low frequency gravity wave background. 




Log[GW Frequency (y ')] 

Fig. 4. — Simulated ULF gravity wave backgrounds from fsi = 1 mergers and loss cone limited 
stellar encounter (MV) mergers. Bright nearby sources rise above the background. Approximate 
sensitivities from anticipated pulsar timing experiments are shown as dashed lines (see text). 

If one makes Ngf^g ~ 20 arrival time measurements per year over a period Tg^g ~ lOy with an 
accuracy of St ~ IS-q/is, then one can place a limit on the strain amplitude of a passing gravity 
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wave with period Py, of roughly 

h ~ |^(iVo6.r„6.)-^/2 ~ 2 X 10-i55_6(iV2oTio)-i/2p-i. (21) 

For this estimate to hold we must have T^hg ^ Pyj, so that fitting pulsar parameters does not 
significantly absorb any gravity wave signal (Blandford et al. 1984). This also assumes that 
the pulsar timing residuals are distributed as white noise. According to this estimate, intensive 
long-term timing programs can reach the sensitivity needed to detect the brightest BHB gravity 
wave sources at periods near lOy. 

Present timing results (KTR) show that PSR1937-I-21 which has been monitored for over eight 
years shows significant unmodeled timing noise, presumably due to rotational variations intrinsic 
to the pulsar. PSR1855-I-09, on the other hand, shows random variations at its ~ 0.8//S arrival time 
accuracy with over seven years of timing; these variations are consistent with perturbations arising 
from instability of the best terrestrial clocks. PSR1855+09 residuals are constraining gravity wave 
sources at an amplitude of /t 5 x 10~^^ (Fig- 4). Thus with the best present data, we estimate a 
chance ~ 0.001 of detecting a merging BHB. However, with a factor of ~ 5 increase in sensitivity 
(and a slightly increased experiment duration) we can anticipate detecting the brightest BHB 
gravity wave sources, for reasonable population assumptions. Whether this sensitivity increase can 
be effected is uncertain. In the last few years, timing programs have been initiated on several new 
pulsars {e.g. PSR J1713+0747, PSR J0437-4715) that provide sub-|Lts timing residuals. However, 
to take advantage of this precision it will be necessary to obtain improved atomic clock standards 
or time one pulsar against another. Finally, to reach interesting detection sensitivities we require 
the arrival time residuals to integrate down as white noise over the ~ lOy periods. Clearly, intrinsic 
instabilities in PSR1937-I-21 prevent this; fortunately timing noise appears to correlate with period 
derivative and several of the new pulsars have period derivatives 10 - 20 times smaller than that 
of PSR 1937+21. It should be noted that the required single order of magnitude improvement 
represents a much better prospect than most other gravitational radiation search techniques! 

In summary, we have developed a model of the inspiral of two massive black holes in a galaxy 
core driven by stellar encounters. We recover some of the behavior described by earlier workers, 
though the large eccentricity needed to circumvent the loss cone catastrophe in our calculations 
must be present ab initio. Our simplified sum allows computation of this process in a range of 
cluster cores, and we see that simple stellar encounter dissipation is only effective in l%-3% of all 
BHB-producing merger events. Nonetheless, it seems likely that other processes will ensure that 
inspiral occurs in the majority of cases. Turning to the rate of galaxy mergers, and the fraction 
of merging cores containing high mass black holes, we estimate the event frequency and typical 
gravity wave amplitude expected from a cosmological population of merging cores with central 
BHBs. Simulating the gravity wave spectrum produced by this population shows that, while no 
detection of ULF GR sources is expected to date, moderate improvements in present sensitivities 
will make detection of waves from binary sources with periods of ~ lOyr possible via timing 
observations of millisecond pulsars. Detection of correlated gravity wave signals in the arrival 
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times of several pulsars would, of course, constitute an exciting confirmation of this astrophysical 
class of gravity wave sources. However, even upper limits a factor ~ 5 lower than present bounds 
can constrain the population and merging behavior of massive BHBs throughout the universe. 

It is a pleasure to acknowledge useful discussions with Andrew Gould and Douglas Richstone, 
and helpful comments by the referee. RWR was supported in part by an Alfred P. Sloan fellowship, 
and MR in part by a fellowship from the National Science and Engineering Research Council of 
Canada. 
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